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Broadband light sources play essential roles in diverse fields, such as high-capacity optical communications, 
optical coherence tomography, optical spectroscopy, and spectrograph calibration. Although a nonclassical 
state from spontaneous parametric down-conversion may serve as a quantum counterpart, its detection and 
characterization have been a challenging task. Here we demonstrate the direct detection of photon numbers 
of an ultrabroadband (110 nm FWHM) squeezed state in the telecom band centred at 1535 nm wavelength, 
using a superconducting transition-edge sensor. The observed photon-number distributions violate 
Klyshko's criterion for the nonclassicality. From the observed photon-number distribution, we evaluate the 
second- and third-order correlation functions, and characterize a multimode structure, which implies that 
several tens of orthonormal modes of squeezing exist in the single optical pulse. Our results and techniques 
open up a new possibility to generate and characterize frequency-multiplexed nonclassical light sources for 
quantum info-communications technology. 



Discriminating the number of incoming photons shot-by-shot plays an indispensable role in quantum 
information processing (QIP) 1 3 , and quantum communication 4 6 . Photon number resolving detection 
(PNRD) is a highly nonlinear process, and hence can be used, combined with nonclassical light sources, 
for quantum state preparation 7-8 and to implement quantum photonic gates 1-3 . Devices for PNRD, such as 
avalanche photodiode (APD) for the time-multiplexing scheme 9 and superconducting transition-edge sensor 
(TES) 10 ~ 13 , have broadband sensitivity, covering most of the spectra of spontaneous parametric down-conversion 
(SPDC), which is an established source of nonclassical light 1415 . SPDC sources inherently generate broadband 
multimode nonclassical states 16-18 . If appropriate mode engineering can be made in those multimode states, the 
implementation of QIP 19,20 can potentially be more scalable. 

Multimode nonclassical light sources and their characterization have been studied extensively these years. In 
the frequency domain, quantum continuous variable (CV) correlations in the frequency comb were generated 
using a single optical parametric oscillator (OPO) and characterized for many sideband modes by homodyne 
detectors with designated local oscillator (LO) modes 21,22 . In the time domain, CV entanglement over massive 
temporal modes, was created by using two OPOs and a delayed interferometer 23 . In the spatial domain, CV 
entanglement over multiple transverse modes was generated from an OPO and characterized by homodyning 24-25 . 

Besides these studies with CVs, discrete variable (DV) nature, i.e., photon-number statistics, of multimode 
quantum states have been investigated and characterized in this decade, by using PNRDs. The direct observation 
of nonclassical photon-number oscillation from a type-I SPDC source was demonstrated using a visible light 
photon counter (VLPC) 26,27 . The higher order photon-number correlations between the signal and idler beams of 
a twin beam from a type-II SPDC source were observed, and the influence of multimode structure to photon- 
number statistics was analysed, using the APD-based time-multiplexing PNRD 28 ' 29 . Extending this idea, recon- 
struction of the mode weight structure from the photon-number correlations was also demonstrated for the states 
containing three modes 30 . All of these works with CVs and DVs have been performed at the visible or near- 
infrared wavelengths. 
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Figure 1 | Experimental setup. Schematic diagram of the experimental setup. The fundamental light is guided into a polarization-maintaining fibre 
(PMF, used as a spatial-mode filter) via a fibre-coupler (FC), passes through variable attenuator (Var. att.) consisting of a half wave plate and a polarizing 
beam splitter, and is then injected into a PPKTP crystal for SHG. The SHG output pumps another PPKTP crystal for squeezing. DMls (DM2s) are 
dichroic mirrors acting as a high reflector for the fundamental (pump) light and as a high transmitter for the pump (fundamental) light. SPDC output is 
coupled into an optical fibre and guided to a TES or a spectrometer. 



In this paper, we demonstrate the direct detection of nonclassical 
photon statistics of ultrabroadband squeezed states at the telecom 
wavelength. The squeezed states are generated from a SPDC source 
and their broadband multimode photon-number statistics is mea- 
sured by a highly efficient PNRD based on a titanium TES 11,12 . The 
centre wavelength of the generated state is 1535 nm and the mea- 
sured bandwidth is as broad as 1 10 nm (13.4 THz), covering the S-, 
C-, and L-bands which are major bands for optical fibre commun- 
ication. We characterize the nonclassical property by the following 
two observations: first, we directly measure the photon-number dis- 
tribution of a whole multimode state and confirm its nonclassicality 
by Klyshko's criterion 26-31 . Second, we apply the method proposed by 
Christ et a\} % and reconstruct the mode weight distribution, i.e. the 
degrees of squeezing in each decomposed mode from the higher 
order photon-number correlations, which implies that several tens 
of independent squeezers are contained in our source. The results 
and techniques pave a new simple way of characterizing ultrabroad- 
band nonclassical states and we anticipate it will accelerate further 
research and development on quantum information technologies for 
the telecom fibre infrastructures. 

Results 

Broadband squeezed light source. Our light source is a single-beam 
squeezer in a collinear configuration where the pump, signal, and 
idler beams are in the same spatiotemporal and polarization mode. 
The generation process is mathematically described by the following 
unitary transformation 1832 



S = exp 



A I doj s I dWjf(co s , coi)a' (co s )a' (a>i) +h.c. 



(1) 



Here A denotes the overall efficiency of the squeezing,/(co s , c/jj) is the 
joint-spectral distribution (JSD) in terms of the signal and idler 

angular frequencies, co s and C0i, and flJ(co s ) ^'(cu;)^ represents 

the creation operator of the signal (idler) field in the continuous 
spectrum of co s (fflj). These continuous modes are, however, not 
necessarily a convenient representation to deal with the multimode 
characteristics, because they cannot be decoupled from each other 
for most of practical JSDs. For a JSD which distributes in an 



effectively finite range, one can find a more convenient discrete set 
of decoupled modes 3 " 6 . 

In fact, when the JSD is engineered to be symmetric, which is the 
case here, it can then be decomposed into the following form 



A*/*(o 8 , COi) = ^2 Tk^tipMkiPi)- 



(2) 



in terms of a complete orthonormal set {^(co)} 18 . Note that each 
(j)k(to) represents the shape of the frequency distribution for signal 
and idler. Here is a modal amplitude corresponding to a squeezing 
parameter for mode k. Introducing the field operators for the broad- 
band basis modes 36 , 



Ak = I doj 



the above transformation can be expressed as 



= exp 



E 



r k [A 



1+2 



(3) 



(4) 



The generated state is a multimode squeezed vacuum state as 

\r) = S\0). (5) 

Thus this state only includes even number of photons. In practice, 
however, inevitable losses in the generation, propagation, and detec- 
tion processes cause the vacuum invasion, making the pure state |r) a 
mixed one, whose statistics has odd number components as well. If 
losses are below a certain level, one can still observe even-rich 
photon-number statistics. Direct observation of this even-odd num- 
ber oscillation and multimode analysis on the distribution of r k are 
our tasks here. 

Experimental setup. Figure 1 shows the experimental setup. The 
fundamental light source is a single-longitudinal-mode, passively 
Q-switched, diode-pumped solid-state laser (Cobolt, Tango), 
operating at 1535 nm with a pulse width of 4 ns and a repetition 
rate of 3.1 kHz. This is used as the fundamental light, and is injected 
into a 10 mm long, periodically poled KTiOP0 4 (PPKTP) crystal for 
second harmonic generation (SHG) at 767.5 nm. This second 
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harmonic light is then used as a pump for squeezing. An unconverted 
fundamental light after SHG is used as a guiding beam for fibre- 
coupling to the TES, but effectively eliminated by dichroic mirrors 
when photon counting experiment is implemented. The pump is 
guided into the second type-0 PPKTP crystal employed as the 
squeezer. The PPKTP crystals for SHG and the squeezer are 
mounted on copper stages and temperature-stabilized. The multi- 
mode squeezed vacuum state is fibre-coupled and guided to the TES 
or to a spectrometer. Our TES is based on a titanium superconductor 
whose size is 10 X 10 /mi 2 as shown in the reference 11,12 , and is set 
inside a dilution refrigerator (TS-3H100-PT, Taiyo Nippon Sanso). 
Output waveform from the TES is amplified by a SQUID device and 
room-temperature electronics (Magnicon). Each waveform was 
sampled using a digital oscilloscope (DPO7104, Tektronix), and 
sent to a host computer for noise filtering and data analysis. 

Figure 2a shows raw waveforms from the TES. Total 1 X 10 7 
waveforms are stored and overlaid here. The gating time for each 
detection event was 10 microseconds. Within this gating period, an 
average photon number of TES's dark counts was measured to be 5 X 
10~ 4 per pulse, which corresponds to 50 counts per second. The 
wavelength dependence of our TES's DE is shown in Fig. 2b (also 
see Methods). It is almost flat from 1512 nm to 1580 nm, and then 
gradually drops from 1580 nm to 1620 nm. The absolute DE is = 
90.9 ± 2.4% at 1535 nm. Unfortunately, we were not able to measure 
DE at wavelengths shorter than 1510 nm nor longer than 1620 nm, 
because of a limited wavelength range of the probe light. 

The red line in Fig. 2c is a spectrum of the fundamental light, 
whose linewidth is —0.05 nm. The green dotted line is a spectrum 
of the SPDC measured with a spectrometer (Acton) for a 30 sec 
accumulation time. Measured full-width at half maximum 
(FWHM) is roughly 130 nm. The blue line is a theoretically derived 
spectrum using the present experimental parameters (see Methods), 
which is symmetric about the fundamental light spectrum. The 
asymmetry seen in the measured spectrum (green dotted line) is 
an artifact caused by the DE fall-off of an InGaAs photo detector 
in the spectrometer. Actually the spectrometer DE gradually drops 
when the wavelength exceeds 1500 nm. Then the DE starts steeply 
falling at around 1580 nm and reaches almost 0% at around 
1650 nm. The blue theoretical spectrum captures the actual SPDC 
spectrum whose FWHM extends over almost 150 nm in the telecom 
window, covering the S-, C-, and L-bands. 

Nonclassical photon statistics of multimode squeezed states. 

Figure 3 shows the experimental results for the multimode squeezed 
states. Through 3a, b, and c, the histograms in red and blue 
correspond to a pump pulse energy of ~1 pj and —10 pj, 
respectively. Figure 3a shows the pulse height distribution of our Ti- 
TES output waveforms. A voltage range of —20 to 60 mV is divided 



into 800 bins. The total number of events was 1 X 10 7 for each data 
set. The vertical axes in the main graphs are numbers of counts shown 
in log scale, while those in the insets are in linear scale. Every peak in 
each distribution was clearly separated, and was fitted with a Gaussian 
function in order to determine thresholds for calculating photon- 
number probabilities. The FWHM of the each Gaussian peak 
corresponds to an energy resolution and was calculated to be about 
0.2 eV, which is smaller enough compared with a photon energy 
(—0.8 eV) of the monochromatic fundamental light at 1535 nm. 

Figure 3b insets show the photon-number distributions obtained 
from the pulse height distributions without any correction of losses. 
In each condition, the results clearly show super-Poissonian distri- 
bution, whose feature could be confirmed by the Fano factors, which 
are greater than the unity, yielding 1.462 ± 0.002 (red) and 1.494 ± 
0.001 (blue), as well as by g* 2) values, 42.746 ± 0.195 (red) and 4.978 
± 0.008 (blue). Then, the overall DE including the fibre-coupling 
efficiency can be evaluated by comparing the single- and two-photon 

2^2/^1 

probabilities, according to the formula, ;; DE = — 27 . The 

l + 2P 2 /"i 

overall DE was thus estimated to be >? DE = 48.4 ± 0.2%. Note the 
dark counts slightly changed P u and we therefore subtracted its 
contribution from the raw statistics for calculating DE. 

The main graphs in Fig. 3b are reconstructed photon-number 
statistics using the iterative maximum-likelihood (ML) estimation 
algorithm 37 . Through the ML estimation process, we used the pos- 
itive operator-valued measure (POVM) of our TES described as 

m M 

follows 38 , n m = exp( — d) — ^2 B(m-i)n\ n )( n \> where 

i= 0 " n — m — i 

B {m - l)n =^ m n _.y i ™ E , (l-T lDE ) n -( m - i) . The POVM elements 

are normalized in order to fulfill S m Yl m = I. TES is thus modelled 
as a linear photon counter which detects m photons from n input 
ones with the finite detection efficiency )7 DE , and with the dark counts 
per pulse d. M is a cutoff photon-number chosen to be large enough 
in order to prevent estimation errors by truncating higher photon- 
numbers. We chose M = 10 here. The reconstructed statistics were 
fully converged after 1 X 10 6 iterations of the ML estimation algo- 
rithm. Note, in our case, the linear inversion of the binomial distri- 
bution B( m -,)A: resulted in unphysical, negative photon probabilities 
because of alternating signs in B^_^ k components 39 . The recon- 
structed two-photon probability P 2 in red (upper panel) is about 
1%. On the other hand, the reconstructed two- and four-photon 
probabilities for the stronger pumping condition (lower panel in 
blue) resulted P 2 ~ 10% and P 4 — 0.6%, respectively. Thus the 
reconstructed probabilities exhibit even-rich photon-number statist- 
ics due to the pairwise photon generation of squeezing. 
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Figure 2 | TES specifications and SPDC spectrum, (a) Typical waveforms from the TES. Total 1 X 10 7 waveforms are overlaid. The vertical direction 
is divided into bins and voltage counts are summed with respect to each bin. Each waveform contains 1000 points (horizontal direction), 
(b) Wavelength dependence of TES's detection efficiency from 1512 nm to 1620 nm. (c) Spectrum of the SPDC. The red line shows spectrum of the 
fundamental laser source. The green dotted and blue lines are measured and theoretical spectra of SPDC, respectively. 
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Figure 3 | Photon statistics of multimode squeezed states, (a) Pulse height distributions of TES output waveforms for two different pump conditions: 
— 1 pj for the upper panel and —10 pj for the lower panel, respectively. The main graphs are in log scale, while the insets are in linear scale, 
(b) Reconstructed photon-number distributions with the loss and dark conts compensated. The insets are raw photon-probability distribution in log 
scale for the vertical axes, (c) Klyshko's figure of merit in log scale. All the error bars in (b) and (c) were statistical errors calculated by linear error- 
propagation method. We assume that photon-counting events obey Poissonian statistics. 



We also evaluated the nonclassicality using Klyshko's criterion 26,31 . 

P P 

Klyshko defined the figure of merit K„ = (n + 1) = 

nP\ 

(n — 1, 2, . . .) and showed that for any {P n } originating from classical 
states, such as coherent state or thermal state, K n s 1 should be 
satisfied for all n 31 . In other words, if this inequality is violated for 
any of K„, the distribution {P„} is necessarily from a nonclassical 
state. In Fig. 3c, K„ were calculated from the raw photon distribution 
(Fig. 3b insets), and the classical limit (green horizontal line) is clearly 
violated (K„ < 1) for even-photon numbers (n = 2, 4, and 6), while 
all the odd photon-numbers do not violate the limit. Ideally, 
squeezed source consists of even-photon numbers and therefore all 
the K„ for odd numbers must go infinite. There were in practice some 
contributions of odd-photon numbers caused by the finite DE, 
resulting in the finite K„ for odd photon-numbers. 

Mode analysis of multimode squeezed states. Using the above 
PNRD data, we evaluate experimental and g* 3 ', and analyse 
multimode structures of our squeezer by applying the method in 
the ref. 18. Figure 4a shows experimental g 0 ^ and g< 3) calculated 
from the photon-number statistics using mth factorial moment 



given by g 



(m). 



Am) 



(«exp)" 



(m 



2 or 3) where 



Am) 



i! 



-P„ and (n exp ^) = nP n A0 . The red and blue 



(n — m). 

n=m v ' «— 0 

solid lines are theoretical plots for the single-mode squeezed state 

given by p' 2 ' = 3 + -j—r , e' 3 ' = 15+-— . Because p*" 1 ' is a loss- 
(n) (n) 

tolerant measure, we use (n) = (n e xp)/'?DE m order to compare the 

experimental g* 2) andg <3) to the theory. As expected from the fact that 

our squeezed states are highly multimode, our experimental g^ 2) and 

g* 3 ' show clear deviations from the single-mode theory. This can be 



contrasted to an experimental result on single-mode squeezed states 
from a PPKTP, which showed that the experimental £ 2) agreed well 
with the single-mode theory 41 . 

It is generally a nontrivial task to carry out the Schmidt decom- 
position analytically. For certain classes of states, however, it can be 
made. A typical class is a two-dimensional (2D) JSD f(a) s , co ; ) in a 
Gaussian form, which is actually the case here. In this case, the eigen 
functions and eigen-value distribution of Eq. (2) are the Hermite 
polynomials and the thermal distribution, respectively, as shown in 
Appendix in the ref. 17. Theng* 2 ' and^ 3) can be expressed as follows 18 
(also see Methods), 



B 2 



(6) 



(7) 



with £ 4 (//) = 



and£ 6 (^() = 



-l+/i 2 



- , where B is an optical 



l + H 2 ° vrv l+/' 2 -ir 

gain coefficient defined by = Ba^ with the normalized mode 

weight At (^2k ~ • ^ ne sm gl e parameter \x determines the ther- 
mal mode distribution as A* = \J 1 — fi 2 /.i k . The slope of 
£mdti- vs -£multi can thus be defined as 

e?Gu)=3 + 6£ 4 Gu). (8) 

The red crosses in Fig. 4b show theg* 3 ' values as a function ofg* 21 . By 
fitting the slope of g^'-vs-g 12 ' to S(fl) using the linear least squares 

fitting technique, we derived s(^/i exp ^ =3.241 + 0.151. This corre- 
sponds to ^ exp = 0.961 ± 0.024. Once ^( exp is determined, each £ exp is 
obtained by eqs. (6) and (7) using corresponding g (2) andg* 3 ', yielding 
0.504 (for ~ 10 pj) and 0.157 (for ~1 pj). Note that £ exp can be 
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independently determined by the measured average photon-number 
via the relation B exp = yj («exp) /^de which implies 0.506 (for ~ 
10 pj) and 0.148 (for ~ 1 pj) agreeing well with the above results. 
Thus we can finally obtain the mode distribution of the squeezing 
parameter (=B e xp/-A:)- The amount of squeezing for each can be 
calculated via — 10 log 10 (exp(— 2r k )) [dB]. Reconstructed squeezing- 
level distributions (k = 0, 1, .. ., 40) are shown in Fig. 4c for a pump 
energy of —10 pj. The effective mode-number K of the squeezed 
vacua contained in the reconstructed distribution can be calculated 

using the relation K= 1 /£ 4 (/j), yielding K^ exp j = 25^*}. 

It is worth to compare /.i exp with its theoretical predictions. From 
the JSD characterized by the phase-matching condition of our non- 
linear crystal and an envelope of the pump pulse, we can calculate the 
theoretical f.i as |ij SD ~ 0.99844 (see Methods for detailed calcuation) 
which is slightly deviated from /( exp . We also estimated the origin of 
the uncertainty of /( exp by the Monte Carlo simulations with 1 X 10 7 
trials (same as that of the experiment. See Methods for details of the 
simulation). Figure 4d compares the experimental and simulated 
g^'-vs-g* 2 ' plots, showing that those agree very well. Then from the 
simulated result, we obtained cS(/( sim ) by linear least squares fitting 
which yieldsS(^ sim ) =3.0164 + 0.1288 and /.i sim = 0.9973 ± 0.0212. 
We conclude that the uncertainty of our experimental data is mainly 
due to the statistical errors. For further investigation of the discrepancy 



between /( exp and /Xjsd> the setup should be further stabilized to allow a 
longer data acquisition for reducing the statistical errors. 

Discussion 

In this work the nonclassical photon-number distributions as the 
even-odd photon-number oscillations were observed in an extremely 
broadband of 1 10 nm FWHM. This bandwidth is, to our knowledge, 
the broadest one whose nonclassicality was directly detected. The 
bands are technologically important ones, namely the telecom S-, 
C-, and L-bands. The mode reconstruction technique theoretically 
developed in the ref. 18 was applied, which inferred that the several 
tens of orthonormal modes were contained in every single optical 
pulse. 

There are interesting open issues along with this work. First, the 
method used here is just for evaluating the mode weight distributions 
under certain assumptions on the SPDC spectral correlation. More 
sophisticated methods on mode analysis using direct detection with 
PNRDs are desired to be developed, as can be implemented by using 
homodyne detection 21,22 . Second, low-loss mode-separating devices 
should also be developed, which may include passive low-loss band- 
pass filters or recent proposals of quantum frequency conversion 42,43 , 
and combined with highly efficient PNRDs to detect and analyse 
mode characteristics of nonclassical states. Third, when a shaped 
broadband LO pulse 44 is introduced in front of a PNRD, possibly 
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combined with mode-separating devices, one can implement a dis- 
placement operation before a PNRD which could realize a phase- 
sensitive PNRD. A displaced squeezed state is known to show unique 
phase-dependent oscillation in photon-number distributions 45,46 , 
which has never been demonstrated yet. It also provides a means 
to implement quantum receivers to overcome the conventional lim- 
its of optical communications 4-6,47 . 

Finally broadband light sources play essential roles in diverse fields 
in the classical domain, such as high-capacity optical communica- 
tions 48 , optical coherence tomography 49 , optical spectroscopy 50 , and 
spectrograph calibration 51 . Supercontinuum sources in the vicinity of 
the 1550 nm telecommunication wavelength 52 are of particular 
importance for wavelength division multiplexing in fibre-optic com- 
munications 53 . Our ultrabroadband source and detection scheme 
shown here could also provide a way to utilize the multimode non- 
classical states for realizing such capabilities in the quantum domain. 

Methods 

Calibration of TES's detection efficiency depending on the wavelength. We used a 
continuous-wave (CW), wavelength-tunable ECDL (TSL-510, Santec) in order to 
measure the DE's wavelength dependence. The CW light beam from the ECDL was 
connected to the fibre-coupler. One of the outputs from the coupler was guided to our 
TES after passing through the variable attenuator. This power to the TES was heavily 
attenuated (—90 dB) such that the average power was set to be below 1 fW, and was 
monitored via the other output port of the coupler using a power meter, whose 
wavelength sensitivity was precisely calibrated. The coupler's splitting ratio and the 
degree of the attenuation were also precisely measured using the same power meter. 
Because photon arrival was random for the CW source, total counts from the TES was 
collected by the threshold detector (SR400, Stanford Research Systems). Photon 
counts from the TES were a few kHz, and the raw waveforms were not overlapped 
each other. We scanned the wavelength at each 3 nm step from 1512 to 1620 nm, and 
obtained the wavelength dependence of our TES's DE as shown in Fig. 2b. 

The second- and third-order correlation functions for the multimode squeezer. In 

this multimode "single-beam squeezer" case, the second- and third-order correlation 
functions are described in the ref. 18 as follows, 



, = 1 + 2 



E t sinh 4 (r t ) 1 
'[£ t sinh 2 (r t )] 2 I>inh 2 ( rt )' 



, = 1 + 6 



Y, k sinh"(r t ) 



E t sinh 6 (nQ 



[E t sinh 2 (r t )] [E t sinh 2 (r t ) 
3 l 6 EtSinh 4 (rQ 

E t sinh 2 (r t ) [E*stah 2 (r0] 3 ' 



(9) 



(10) 



a(aj s , ojO = sinc(LAfc/2) with the crystal length L and phase mismatch Ak = k p — k s — 
ki — 2n/A. A denotes a grating period of the PPKTP crystals and 24.2 /im in our case. 
Assuming that the transversal phase -matching is perfect, the longitudinal one can be 
described as Ak = k p<z — k StZ — k itZ — 2n/A, because we used collinear configuration 
where pump, signal, and idler are all in the same polarization. Note that the beam 
propagating direction is taken to be along z-axis and fc^fco^J — (O^nJ^ co^/c is the wave 
vectors. The frequencies are in a relation of C0j = co p — 0) s = 2a> 0 — cq s . The 
corresponding Sellmeier equation for n z can be found in the ref. 55. The blue line in 
Fig. 2c was theoretically derived using the above formulas and parameters. 

The theoretical Kjsd can be calculated from the JSD using the singular value 
decomposition (SVD) technique 29 . The JSD was obtained from the above parameters 
and the pump linewidth of — 0.05 nm, and then discretized into a 16000 X 16000 
square matrix with a frequency resolution of 2 GHz and with a frequency range of 
179.3 THz(~ 1672 nm)to211.3 THz(~ 1419 nm). We carried out SVD for the JSD 
matrix using MATLAB, and obtained Kj SD 640. The corresponding \i is ^j SD 
0.99844. 

Theoretical photon probability of multimode squeezed states for Monte Carlo 
simulations. In order to analyse the discrepancy between the experimentally 
obtained /t exp and the theoretically obtained /(jsd' we carried out Monte Carlo 
simulations to illustrate how finite data size affects S(/()-fitting errors. Before the 
simulations, we first calculate theoretical photon probabilities for the multimode 
squeezed vacua. The density matrix of the multimode vacuum state should be 
expressed as a tensor product of the decomposed squeezed vacuum states, whose 
squeezing parameters are characterized by — B exp A k ) . From this density matrix, one 
can derive the M-photon generation probabilities of the source. Though n is 
distributed to infinity, we truncated it at six photons which is a reasonable 
approximation to simulate our experimental situation in Fig 3b. The n-photon 
generation probabilities up to n = 6 are given by 



1 K ' 

ptot _ _ rr P 1 

A/ A: = 0 



,(*) 



1 K 

= — Vp co 



j K K 



j K! K' K' 



k=o m=fc+iR=m+l 



(16) 



where '=2/(2 + ^), Pf ) = r 2 k /{2 + rQ, and Af= ^ p T( n = 0,2,4,6). Here we 
assume that the r k is small enough and each mode only emits vacuum or two-photons. 
This assumption is valid when the mode number K' is large enough. For example, in 
the case of four-photon emission, the number of cases for two-photon generations 
from the two different modes (k-C 2 ~ 0{K' 2 )) are much greater than those for four- 
photon emission from only one mode GrQ ~~ O(K')). Then the Monte Carlo 
simulation was performed along with eq. (16). 



When a JSD takes a 2D Gaussian shape, A k forms a thermal distribution defined by 
a single parameter fi, as Xk — \Jl — / 1 2 \t . When the r k is small enough (r^-Cl), that is 
valid in our experimental conditions, the following approximations are valid as 

1 1 11 

£ZLo sinh 2 (r k ) ~ ULy k ~ W " * ' 



(11) 



Z^fc = 0 r k 



[Er=oSinh 2 (r,)] 2 E 



:£> 4 = £ 4 ( A ,), 



(12) 



E t =o sum (r k ) Et-o h _ \- ,6 _ r , lA 
[Ei:_o smh [Et=0 r tJ t=0 

Then, eqs. (9) and (10) can be simplified into the eqs. (6) and (7)*g*fA ti an dgmiiiti can 
be exactly calculated from the experimental photon statistics. Furthermore, from eqs. 
(6) and (7), g^L; can be expressed as a function of g^ ti as 



«£w = (3 + 6£ 4 ( A1 )ki 2 l -2-6£ 4 ( /i ) - 12£ 4 (rt 2 + 8£ 6 ( M ). 
The slope of g^-vs-g^ is finally defined as eq. (8). 



(14) 



Specification of PPKTP crystal and SPDC spectrum. The JSD can be a product of 
two quantities as 

f(co s ,cOi) = ot(co s ,a>i)^(f(} s ,cOi), (15) 

where a(oj s , describes pump -frequency envelop, and <j>{co s , oj f ) describes 
phase -matching function for the SPDC in a nonlinear crystal 54 . Explicitly 
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